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We study, through large scale stochastic simulations 
using the noise reduction technique, surface growth via 
vapor deposition e.g. molecular beam epitaxy (MBE), 
for simple nonequilibrium limited mobility solid-on-solid 
growth models, such as the Family (F) model, the Das 
Sarma-Tamborenea (DT) model, the Wolf- Villain (WV) 
model, the Larger Curvature (LC) model, and other re- 
lated models. We find that d=2+l dimensional surface 
growth in several noise reduced models (most notably the 
WV and the LC model) exhibits spectacular quasi-regular 
mound formation with slope selection in their dynamical 
surface morphology in contrast to the standard statisti- 
cally scale invariant kinetically rough surface growth ex- 
pected (and earlier reported in the literature) for such 
growth models. The mounding instability in these epitax- 
ial growth models does not involve the Ehrlich-Schwoebel 
step edge diffusion barrier. The mounded morphology in 
these growth models arises from the interplay between the 
line tension along step edges in the plane parallel to the 
average surface and the suppression of noise and island 
nucleation. The line tension tends to stabilize some of the 
step orientations that coincide with in-plane high sym- 
metry crystalline directions, and thus the mounds that 
are formed assume quasi-regular structures. The noise 
reduction technique developed originally for Eden type 
models can be used to control the stochastic noise and 
enhance diffusion along the step edge, which ultimately 
leads to the formation of quasi-regular mounds during 
growth. We show that by increasing the diffusion sur- 
face length together with supression of nucleation and de- 
position noise, one can obtain a self-organization of the 
pyramids in quasi-regular patterns. The mounding insta- 
bility in these simple epitaxial growth models is closely 
related to the cluster-edge diffusion (as opposed to step 
edge barrier) driven mounding in MBE growth, which has 
been recently discussed in the literature. The epitaxial 
mound formation studied here is a kinetic-topological in- 
stability (which can happen only in d=2+l dimensional, 
or higher dimensional, growth, but not in d=l-|-l dimen- 
sional growth because no cluster diffusion around a closed 
surface loop is possible in "one dimensional" surfaces), 
which is likely to be quite generic in real MBE-type sur- 
face growth. Our extensive numerical simulations produce 
mounded (and slope-selected) surface growth morpholo- 
gies which are strikingly visually similar to many recently 
reported experimental MBE growth morphologies. 



I. INTRODUCTION 

Crystal growth, particularly high-quality epitaxial thin 
film growth, is one of the most fundamental processes im- 
pacting today's technology 0. A major issue in interface 
growth experiments is to have continuous dynamical con- 
trol over the deposition process, such that interfaces with 
certain desired patterns can be obtained. For example, 
while in thin film epitaxy it is desirable to obtain smooth 
surfaces, in nanotechnology it is also important to be 
able to create regular, nanoscale structures with well de- 
fined geometry such as quantum dots, quantum wires, 
etc. Growth is usually achieved by vapor deposition of 
atoms from a molecular beam (molecular beam epitaxy, 
or MBE). Thus, in order to be able to design a controlled 
deposition process, it is of crucial importance to under- 
stand all the instability types (which destroy controlla- 
bility) that may occur during growth. In the present 
paper we concentrate on MBE growth. There are sev- 
eral types of instabilities in MBE among which we men- 
tion the Ehrlich-Schwoebel (ES) 0] instability which is of 
kinetic-energetic type. Ehrlich-Schwoebel (ES) barriers 
in the lattice potential induce an instability by hindering 
step-edge atoms on upper terraces from going down to 
lower terraces Q which in turn can generate mounded 
structures during growth. The ES instability is thought 
to be an ubiquitous phenomenon in real surface deposi- 
tion processes, and as such, is widely considered to be 
the only mechanism for formation of mounds in surface 
growth. Another kinetic-energetic instability, which does 
not involve any explicit ES barriers, was discussed by 
Amar and Family ||. This instability involves a neg- 
ative barrier at the base of a step, and thus is due to 
a short range attraction between adatoms and islands. 
Short range attractions generically lead to clustering in 
multiparticle systems, a property which in the language 
of MBE translates into formation of mounds. Note that 
this attraction-induced instability and ES instability are 
essentially equivalent, and both could occur in d=l+l 
or 2+1 dimensions. Both the attractive instability and 
the ES barrier instability cause mounding because atoms 
on terraces preferentially collect at up-steps rather than 
down-steps leading to the mounded morphology. In both 
cases, there is no explicit stabilizing mechanism, and the 
mounds should progressively stiffen with time leading 
to mounded growth morphology with no slope selection, 
where the mound slopes continue to increase as growth 
progresses. It is, of course, possible to stop the monotonic 
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slope increase by incorporating some additional mecha- 
nisms (e.g. the so-called "downward funneling" where 
deposited atoms funnel downwards before incorporation) 
which produce a downhill mass current on the surface 
and thereby opposes the uphill current created by the 
ES barrier. There is, however, no intrinsic slope selec- 
tion process built in the ES barrier mechanism itself. 




FIG. 1. Mounded morphology created by SED. This is a 
simulation of the LC model as described in the text. 



The purpose of the current article is to discuss a com- 
pletely different kinetic-topological mechanism, qualita- 
tively distinct from the ES barrier induced mounding, 
which leads to spectacular mound formation (see Fig. |l| 
for a colorful example) in MBE growth without involving 
any ES barriers (or the closely related island-adatom at- 
traction mechanism) whatsoever. Based on the direct 
stochastic numerical simulations of a large number of 
MBE-related nonequilibrium limited mobility solid-on- 
solid epitaxial growth models studied in this paper we 
believe that epitaxial mounding of the type discussed 
here is quite generic, and may actually apply to a va- 
riety of experimental situations where mounded growth 
morphologies (with slope selection) have been observed. 
The mechanism underlying the mounding instability dis- 
cussed in this paper is closely related to that proposed in 
two recent publications JJH although its widespread ap- 
plicability to well-known limited mobility MBE growth 
models, as described here, has not earlier been appre- 
ciated in the literature. The mounding instability dis- 
cussed in this paper, being topological in nature, can hap- 
pen only on physical two-dimensional surfaces (d=2+l 
dimensions growth) or in (unphysical) higher (d>2+l) 
dimensional systems, in contrast to ES type instabilities 
which, being of energetic origin, are allowed in all dimen- 
sions. 

The instability most recently discovered by Pierre- 
Louis, D'Orsogna, and Einstein M, and independently by 



Ramana Murty and Copper |5| is entirely different from 
ES instabilities, because it is not of kinetic-energetic na- 
ture (involving potential barriers) but rather of a kinetic- 
topologic nature. This instability is generated by strong 
diffusion along the edges of monoatomic steps, and its 
net effect in two or higher dimensions is to create quasi- 
regular shaped mounds on the substrate (see Fig. |l|). 
For simplicity of formulation we will call this instability 
the step-edge diffusion, or SED instability. What comes 
as a bit of a surprise is that in spite of its simplicity 
as a mechanism, it has not apparently been recognized 
that the SED mechanism may actually be playing a domi- 
nant role in many epitaxial mounding instabilities. There 
are two reasons for this, one is probably because the ES 
instability has inadvertently and tacitly been accepted 
as the only mechanism for mounding, and therefore the 
mounded surface morphologies coming from experiments 
have been apriori analyzed with the ES mechanism in 
mind, and the second reason is that SED is not easy to 
see in simulations unless other conditions (related to re- 
duced noise) are met in addition to the existence of edge- 
current, which we will discuss in details in the present 
paper. 

In this paper we discuss the SED instability, and an- 
alyze its effects on various growth models, giving both 
a simple continuum and discrete description. The SED 
is shown to exist in a large number of discrete, limited 
mobility growth models, and we identify the conditions 
under which SED-induced epitaxial mounding is mani- 
fested. We study the instability by producing dynamical 
growth morphologies through direct numerical simula- 
tions and by analyzing the measured height-height cor- 
relation functions ](| of the simulated morphologies. This 
paper is organized as follows: in Section II we define a 
number of limited mobility growth models (Subsection 
II. A) studied here, we briefly review the noise reduc- 
tion technique (Subsection II. B), then present a series 
of simulated growth morphologies, both from early and 
late times, with and without noise reduction (Subsection 
II. C); in Section III we describe and discuss the SED 
mechanism, by first presenting a continuum description 
in terms of local currents (Subsection III. A) and its dis- 
crete counterpart with particular emphasis on the Wolf- 
Villain and Das Sarma - Tamborenea models (Subsec- 
tion III.B); in Subsection III.C we introduce the notion 
of conditional site occupation rates to analyze the effects 
of noise reduction; Section IV is devoted to the prop- 
erties of the mounded morphologies with the emphasis 
on the relevance of the height-height correlation function 
(Subsection IV. A), where we note possible comparison 
with experiments and discuss the relevance of the global 
diffusion current in the simulations for the mounding in- 
stability. 
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II. MORPHOLOGIES OF DISCRETE LIMITED 
MOBILITY GROWTH MODELS 

A. Models 

All the models studied here are dynamical limited mo- 
bility growth models in which an adatom is allowed to dif- 
fuse (according to specific sets of diffusion rules for each 
model) within a finite diffusion length of I sites. The dif- 
fusion process is instantaneous and once the adatom has 
found its final site, it is incorporated permanently into 
the substrate and can no longer move. These models are 
also based on the solid-on-solid (SOS) constraint where 
bulk vacancies and overhangs are not allowed. Desorp- 
tion from the growth front is neglected. 




LC model 



FIG. 2. Schematic configurations defining the diffusion 
growth rules in one dimension for various dynamical growth 
models discussed in this paper. Here 1 = 1. 

With this SOS constraint, growth can be described on 
a coarse-grained level by the continuity equation 

dh 

^ + v-j = ra + 7 ? (x,i), (i) 

where /i(x, t) is the surface height measured in the growth 
direction, j is the local (coarse-grained) surface diffusion 
current, F is the incoming particle flux, Vl — a±a^ is the 
surface cell volume with aj_ and a|| being the lattice con- 
stants in the growth direction and in the ci-dimensional 
substrate, respectively. The shot noise is Gaussian un- 
corrected white noise with correlator: 

(r/(x, t) V (x',t')) = D5(x - x')6(t - t') (2) 



where the amplitude of the fluctuations D is directly pro- 
portional to the average incoming particle flux through 
the relation: D = Ftt 2 . 

The most general form of this continuity equation Q 
which preserves all the symmetries of the problem can be 
written as: 

dh 

— = v ^ 2 h - A 4 V 4 /i + A 22 V 2 (V/i) 2 + ... + rj, 

dt (3) 

where h from now on indicates the height fluctuation 
around the average interface height (h). Note that the 
unit of time in our simulations is defined through (h); 
we arbitrarily choose the deposition rate of 1 monolayer 
(ML) per second, and all simulations are done with peri- 
odic boundary conditions along the substrate. 

1. Family model 

This is a model introduced by Family (F) Q with 
the simplest diffusion rule: adatoms diffuse to the local 
height minima sites within range I (see the one dimen- 
sional version in Fig. 2). This is a well studied model and 
it is described on a coarse-grained level by the continuum 
Edwards- Wilkinson equation ^ = vV 2 h + r/ and 
hence the model asymptotically belongs to the Edward- 
Wilkinson (EW) universality class. 

2. Das Sarma-Tamborenea model 

The Das Sarma-Tamborenea (DT) model was pro- 
posed ]Tl[ | as a simple limited mobility nonequilibrium 
model for MBE growth. In this model, an adatom is 
randomly dropped on an initially flat substrate. If the 
adatom already has a coordination number of two or 
more (i.e. one from the neighbor beneath it to satisfy 
the SOS condition, and at least one lateral neighbor) at 
the original deposition site, it is incorporated at that site. 
If the adatom does not have a lateral nearest neighbor at 
the deposition site, it is allowed to search, within a finite 
diffusion length, and diffuse to a final site with higher 
coordination number than at the original deposition site. 
In the case where there is no neighbor with higher coordi- 
nation number, the adatom will remain at the deposition 
site. At the end of the diffusion process, the adatom 
becomes part of the substrate and the next adatom is 
deposited. 

It is important to note that adatoms in the DT model 
search for final sites with higher coordination numbers 
compared to the deposition site. The final sites are not 
necessarily the local sites with maximum coordination 
numbers. In other words, in the DT model adatoms try 
to increase, but not necessarily maximize, the local co- 
ordination number. Also, deposited adatoms with more 
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than one nearest-neighbor bond do not move at all in the 
DT model. 



combination of h and its various differential forms, V 2 /i, 
(V/i) 2 , ... etc., where certain combinations are ruled out 
by symmetry constraints: 
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FIG. 3. Morphologies from F model with n r = 1 (top) 
and n r = 5 (bottom) from substrates of size L = 100 x 100 
after 10 6 ML. 

3. Wolf- Villain model 

The Wolf- Villain (WV) model U) is very similar to 
the DT model in the sense that diffusion is controlled 
by the local coordination number. There are, however, 
two important distinctions between the two models. For 
one thing, in the WV model, an adatom with lateral 
nearest neighbors can still diffuse if it can find a final 
site with higher coordination number. The other differ- 
ence is that, adatoms in the WV model try to maximize 
the coordination number, and not just to increase it as 
in the DT model (see Fig. ||). The two models have 
very similar behavior in one dimensional (d=l+l) sim- 
ulations, and there have been substantial confusion re- 
garding these two models in the literature. But in two 
dimensional (d=2+l) systems, after an initial crossover 
period, the two models behave very differently, as will be 
shown in this paper. 

4- Kim - Das Sarma class of conserved growth models 

A number of discrete growth models can be defined 
by applying the approach originally developed by Kim 
and Das Sarma |l3| ] for obtaining precise discrete coun- 
terparts of continuum conserved growth equations. In 
this method, the surface current j is expressed as the 
gradient of a scalar field K, which in turn is written as a 



j(x,t) = -V*T, (4) 

where 

K = u 2 h- X A V 2 h + \ 22 {Vh) 2 + ... (5) 




FIG. 4. Morphologies from DT model with n r = 1 (top) 
and n r = 5 (bottom) from substrates of size L — 100 x 100 
after 10 6 ML. 

Using the above expression for K, the general con- 
served growth equation, Eq. (||), is obtained. Due to the 
conserved nature of growth, one concludes that the local 
particle transport happens in the direction of the largest 
variation of the scalar —K, i.e., in the direction where 
K has the largest decay. This allows for the definition of 
the following Kim-Das Sarma discrete atomistic rules: 1) 
a site i is chosen at random, 2) the scalar K is computed 
on the lattice at site i and for all its nearest neighbors, 3) 
then a particle is added to the site that has the smallest 
value of K; if site i is among the sites with the smallest 
K, then the particle is deposited at site i, otherwise one 
picks a final site with equal probability among those sites 
that have the common smallest K. 
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FIG. 5. Morphologies from WV model with n r = 1 (top) 
and n r = 5 (bottom) from substrates of size L = 500 x 500 
after 10 4 ML. 

For example, the F model (where adatoms relax to the 
sites with local height minima) is described by the linear 
second order growth equation, with K ~ h in Eq. 
Thus the well-known F model obeying the linear EW 
equation is the simplest example of the Kim-Das Sarma 
class of discrete growth models. A well-known member 
of the Kim - Das Sarma class is the Larger Curvature 
model (LC) ]l3| , p^ | where K ~ —V 2 h (which generates 
the linear fourth order growth equation). Adatoms in 
this model diffuse to sites with minimum — V 2 /i, i.e. with 
maximum curvature, hence the name of the model. 

In this paper we also study the nonlinear fourth order 
equation (|jl5|] where: 

K = -X i V 2 h + \ 22 (Vh) 2 (6) 

with A22 being a small constant controlling the instability 
in the equation. Recently, a nonlinear fourth order equa- 
tion with controlled higher order nonlinearities [ p~6|Jl~7| 
has been proposed as a continuum description for the 
DT model. The scalar K in this equation is 

K = -X i V 2 h+^[l-e- c ^ 2 }. (7) 

The last variant we analyze is the the linear sixth order 
continuum growth equation, with K ~ — V h. 
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FIG. 6. Morphologies from DT model in the presence of 
an ES barrier with n r = 1 (top) and n r = 5 (bottom) from 
substrates of size L = 100 x 100 after 10 6 ML. 

B. The noise reduction technique 

In this section we briefly discuss a well-known theoreti- 
cal method used to reduce stochastic noise in the stochas- 
tic growth simulations of surface morphologies, in partic- 
ular its short wavelength component. The short wave- 
length noise component generates an intrinsic width (wi) 
in the surface roughness giving rise to strong nontriv- 
ial corrections to scaling |lS[ ], and typically has different 
scaling properties than the long wavelength quantities, 
such as the dynamical width (w ~ t@) and the correlation 
length parallel to the surface (£ ~ t 1 ^ ). Reducing the in- 
trinsic width through the noise reduction technique, the 
scaling behavior improves dramatically, and the asymp- 
totic critical properties emerge in the simulations [[l8] . 
In the simplest growth models of MBE, two main compet- 
ing kinetic processes are kept and simulated: (1) depo- 
sition from an atomic beam, and (2) surface diffusion of 
the freshly landed atoms until they incorporate into the 
surface according to some diffusion rules. The shot noise 
in the beam is simulated in (1) by randomly selecting the 
landing site, and the surface diffusion performed by the 
freshly landed atoms is simulated in (2) essentially by a 
random walk of I steps biased by the local coordination 
(filled and unfilled sites) landscape of the surface. 
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FIG. 7. Morphologies from WV model in the presence of 
an ES barrier with n r — 1 (top) and n r = 5 (bottom) from 
substrates of size L — 100 x 100 after 10 6 ML. 

The noise reduction method, originally introduced 
p9| , p0| to study the Eden model, proceeds as follows: a 
site is chosen randomly (with uniform probability L~ d ), 
then the /-step random walk ("diffusion") is performed 
according to the local diffusion rules of the model. As- 
sume that the site where the atom would be incorporated 
is site i. Instead of actually depositing the atom at site i, 
first a counter Cj associated with site i is increased. The 
deposition attempt becomes a "true" deposition process 
only if the counter at the attempt site reaches a prede- 
termined threshold value n r . After a true deposition the 
counter is reset to zero, and the process is repeated. The 
integer n r is the noise reduction factor and it quanti- 
fies the level of noise suppression, larger values meaning 
a stronger supression. Obviously n r — 1 simply means 
that there is no noise reduction. The noise reduction 
technique has been successfully used in limited mobility 
MBE growth models |L8|,^l). The effects of noise reduc- 
tion can also be understood as a coarsening process, as 
shown in Ref. Q: noise reduction performs an effective 
coarsening of the shot noise within regions of linear size 
£, which is defined as the linear distance within which no 
nucleation of new islands can form on the terraces. 

An important observation is that when the noise re- 
duction technique (NRT) is used, the 'material parame- 



ters' such as the noise strength, diffusion constant, step 
stiffness, etc. in general become dependent on the level 
of noise reduction M, which is quantified by the noise 
reduction parameter n r . 




FIG. 8. Morphologies from LC model with n r — 1 (top) 
and n r = 8 (bottom) from substrates of size L — 200 x 200 
after 10 3 ML. 

Brendel, Kallabis and Wolf have recently identified 
the dependence of the material parameters on the level 
of noise reduction for kinetic roughening described on 
coarse-grained scale by the KPZ equation, ^ = vV 2 h + 
X(Wh) 2 + Ffl + r\ in the context of layer-by-layer growth. 
In this case it turns out that the parameters v and A 
assume power law behavior: v ~ (n r ) e " and A ~ (n r ) ex , 
where the exponents are not universal and may depend 
on the details of the way in which the model is simulated. 
A simple argument by Kertesz and Wolf shows that 
the noise strength D is weakened by NRT according to: 

D(n r ) oc l/n r . (8) 

Two of the current authors have earlier studied |l^] the 
effect of noise reduction in the growth simulations of lim- 
ited mobility DT, WV, and F models, finding that NRT 
is extremely effective in suppressing crossover and cor- 
rection to scaling problems in these growth models, lead- 
ing effectively to the asymptotic critical exponents rather 
compellingly. An important feature of the current arti- 
cle is to obtain simulated 2+1-dimensional growth mor- 
phologies of the various growth models under NRT. We 
find that NRT is extremely effective in suppressing noise 
and giving rise to the epitaxial mounded morphologies in 
models where SED mechanism is operational. 
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FIG. 9. Morphologies from nonlinear fourth order model 
with n r = 1 (top) and n r = 5 (bottom) from substrates of size 
L = 100 x 100 after 10 3 ML. The coefficient of the nonlinear 
term is A22 = 0.5. 

Often (almost always) the limited mobility models of 
interest in our paper are far too noisy (because of the 
extreme limited role of diffusion in the models) to exhibit 
epitaxial mounding without effective noise reduction, at 
least in treatable simulation times. 

C. Growth morphologies 

In this section, we show simulated dynamical growth 
morphologies of the discrete SOS models described ear- 
lier, with and without NRT. 

The first model is the F model which is known to 
be described exactly by the linear second order growth 
equation and belong to the well known EW universality 
class. As shown in Fig. F morphology with n r = 1 
(without NRT) is already quite smooth. However, when 
noise is suppressed with n r = 5, the layer- by-layer growth 
process persists as long as 10 6 monolayers, which is our 
longest simulation time for the system. With layer-by- 
layer growth, the interface is extremely smooth and flat 
and this agrees with the EW universality definition where 
the growth {(i) and the roughness (a) exponents in 2+1 
dimensions is zero corresponding to smooth growth. 



FIG. 10. Morphologies from nonlinear fourth order model 
with an infinite series of higher order nonlinear terms with 
n r = 1 (top) and n r — 5 (bottom) from substrates of size 
L = 100 x 100 after 10 3 ML. The coefficients of the nonlinear 
term are A22 = 5.0 and C = 0.085. 

The DT model, as shown in Fig. ||, shows kineti- 
cally rough interface in the n r = 1 simulation. But with 
n r = 5, the morphology becomes considerably smoother, 
without any specific mounding pattern. Similar to the 
DT model, the WV model (Fig. |J) with n r = 1 shows 
rough morphology without any special mounding fea- 
tures. But when we suppress the noise, using n r > 1, 
the WV morphology changes drastically and patterned 
mound formations can clearly be seen. This mound for- 
mation in the WV model is intriguing because there is no 
ES barrier involved in the simulation. As it turns out, the 
WV diffusion rule has an intrinsic machanism to create a 
local uphill particle current which is strong under NRT, 
and hence a pyramid-like mounded morphology formed. 
This topic is discussed in the next section. 

We have actually carried out both DT and WV model 
growth simulations |2lJ in the presence of an ES bar- 
rier also. For the sake of comparison with the epitaxial 
mounding in the WV model without any ES barrier (but 
with NRT), as shown in Fig. ||, we present in Figs. || 
and the observed mounded morphologies in the DT 
and WV models in the presence of an ES barrier (both 
with and without NRT). One can clearly see that the ES 
barrier induced mounding instability in Figs. ^ and Q 
is qualitatively different from the epitaxial mounding in 
Fig. ||. While the ES barrier induced DT/WV mounding 
has flat tops, the intrinsic WV mounding of Fig. || is of 
pyramid shape. 
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FIG. 11. Morphologies from linear sixth order model with 
n r — 1 (top) and n r — 5 (bottom) from substrates of size 
L = 100 x 100 after 10 6 ML. 

The linear fourth order model (the LC model) also 
creates mounds without ES barrier in a way similar to 
the WV model. This is shown in Fig. || (and also Fig. 

0) - ...no 

Although the epitaxial mounding in Figs, km and H look 

similar, the mounded morphologies in LC and WV mod- 
els are not identical, however, as the LC is a linear model 
and hence the morphology exhibits up-down symmetry, 
while the WV morphology does not look the same when 
it is turned upside down, and therefore does not possess 
the up-down symmetry. This lack of up-down symmetry 
in the WV model arises from the presence of the nonlin- 
ear V 2 (V/i) 2 term in its dynamics. When we add a small 
instability into the LC equation, by adding the nonlinear 
fourth order term with a small coefficient, the interface 
(shown in Fig. ^) is somewhat mounded even for n r = 1 
simulation. But the surface in this n r = 1 model is so 
rough that a specific mound pattern cannot be discerned. 
When we reduce the noise, using n r = 5, clear mound for- 
mation (without up-down symmetry) can be seen in the 
morphology as one would expect. 

In Fig. [l^ we present our simulation results for the 
foil nonlinear equation (Eq. (Q) with C = 0.085) with 
an infinite number of nonlinear terms. The rationale for 
studying this particular continuum equation with an infi- 
nite number of nonlinear terms of the V 2 "(V/i) 2 ™ form is 
the fact that recently it has been established (l(]|l7| that 
such as infinite order nonlinear dynamical equation is the 



most likely continuum description (at least in d=l+l di- 
mensions, where this issue has so far been studied inten- 
sively) of the discrete DT model. 
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50 100 150 200 
FIG. 12. Contour plots of LC morphologies from 

L = 200 x 200 at 10 3 ML with n r = 8. Lighter (darker) 

shade represents higher (lower) points. Plot (a) corresponds 

to the "higher priority" version while plot (b) is from the 

"equal priority" version. Here 1 = 1. 

Since the DT model does not seem to exhibit (see Fig. 
|j) any epitaxial mounding instability in our simulations, 
it is of considerable interest to see how the infinite order 
nonlinear continuum equation behaves. It is, therefore, 
gratifying to see, as shown in Fig. |l^, that the growth 
morphology in the infinite order nonlinear equation does 
not exhibit any discernible mounding patterns — the 
morphologies depicted in Fig. [l^ look vaguely similar 
to the kinetically rough growth of the DT model shown 
in Fig. ||. It is, however, remarkable that the fourth or- 
der continuum equation (either the linear LC model as 
in Figs. [j] and || or the nonlinear model with a small 
nonlinearity as in Fig. ^|) manifests striking mound for- 
mation, but the infinite order nonlinear equation does not 
manifest any epitaxial mounding. We do not have a pre- 
cise mathematical understanding of this dichotomy other 
than pointing out that the presence of an infinite series 
of nonlinear terms clearly overwhelm the SED instability 
arising from the V 4 /i surface diffusion term. The precise 
understanding of the combined effects of noise reduction 
and growth nonlinearities is not available at the present 
time. Finally, in Fig. [ll] we present the surface mor- 
phologies associated with the linear sixth order growth 
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equation. The mounding pattern can be barely seen in 
the n r = 1 simulation. But the spectacular mound for- 
mation can be seen clearly when we reduce the noise 
(n r — 5). 
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FIG. 13. Contour plots of LC morphologies from 
L = 200 x 200 at 10 3 ML with n r = 8. Lighter (darker) 
shade represents higher (lower) points. Plot (a) corresponds 
to the "higher priority" version while plot (b) is from the 
"equal priority" version. Here I = 2 

These simulated morphologies are very sensitive to the 
local atomistic diffusion rules, as can be seen when we 
make some minor modifications to the LC rule. The hrst 
thing we modify is to make the adatoms more "ener- 
getic". In the original LC model, an adatom searches 
for a site with the largest curvature, and if the original 
random deposition site is one of the sites with the largest 
curvature then the adatom remains at the deposition site. 
So in the original version of the model, an adatom will 
not move if it is already at one of the preferred sites. We 
call this original model the "higher priority" version of 
the LC model, as the deposition sites have higher prior- 
ity than the neighboring sites. In the modified version, 
called "equal priority" version, the adatom chooses its 
final site among all the sites with largest curvature with 
equal probability. So even if the original deposition site 
is one of the largest curvature sites, the adatom may still 
end up being incorporated at a neighboring site if that 
neighbor also has the same curvature as the deposition 
site. The second modification we make is simply to in- 
crease the diffusion length I from nearest-neighbor-only 



diffusion (I = 1) to I > 1 where diffusion rules now extend 
to Z(> 1) sites. 




It is known (Tl]] that increasing diffusion length does 
not change the asymptotic universality class of limited 
mobility growth models as long as the substrate size is 
large enough so that finite size effects remain neglegible. 
But increasing diffusion length can most certainly change 
the early-time morphology of the system by strongly af- 
fecting finite size effects. 

Both versions of LC model show mounded morpholo- 
gies, as can be seen in Fig. |l2[ But it seems that the 
"equal priority" version, the one with more energetic 
adatoms, has mounds with the bases lined in the 110 di- 
rection. In the "higher priority" version the mounds are 
more stable when the bases are in 100 direction. Increas- 
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ing diffusion length in the simulations encourages faster 
coarsening process (Fig. as should be expected, and 
individual mounds in the I — 2 simulations "link" to 
each others more than in the / = 1 simulations. To fur- 
ther investigate the mound coarsening process, we study 
the "equal priority" LC model with I = 2 in long time 
limit and find that the mounds rearrange themselves into 
a regular pattern as shown in Fig. [l4|. This regular re- 
arrangement of the mound morphology, caused by longer 
diffusion lengths as well as the different symmetry prop- 
erties (110 versus 100) of the equal priority and higher 
priority LC models, show the extreme sensitivity of the 
epitaxial mound morphology to the details of the dif- 
fusion process controlling growth — in general, we find 
the SED instability in all our simple models to be quite 
complex in the sense of the mound patterns (and even 
whether mounds exist or not as in WV and DT models) 
seem to be very sensitive to the details of the surface 
diffusion rules. 



III. KINETIC-TOPOLOGICAL "INSTABILITY": 
STEP-EDGE DIFFUSION 

A. A continuum description 

It has been recently recognized || that step-edge dif- 
fusion generates a local uphill surface diffusion current 
during the deposition process, leading to a morphological 
"instability" , i.e. the starting substrate is no longer sta- 
ble after growth. This instability is purely of topologic- 
kinetic nature, and it is not generated by an energetic 
instability such as the Ehrlich-Schwoebel instability. It 
is rather a consequence of the fact that surface diffu- 
sion currents (in two or higher dimensions) are vectorial 
quantities, and while some of their components have a 
smoothening and stabilizing effect in a certain subspace, 
in other directions they may act as instabilities, as is 
shown schematically in Fig. ^|. This figure is a top view 
of a two dimensional substrate on which there are two 
vicinal surfaces. Due to diffusion, the adatoms on the 
flat terraces may reach the edge of the terrace, however 
they are assumed to feel no energetic barriers at the step 
edges. After reaching the edge of the step the adatoms 
will tend to remain attached to it since it is an energeti- 
cally more favorable position (with at least one more ex- 
tra bond when compared to the position on the terrace) . 
The atoms on the step edge may relocate through line- 
diffusion along the edge to energetically even more favor- 
able places (kink sites), thus reducing the differences in 
the local chemical potential. (For simplicity we assume 
no desorption from the growth front.) For example, if 
there is a protuberance on (or a dent in) the step as 
shown in Fig. [l^ (a), the step-edge atoms will preferen- 
tially diffuse towards the base of the hump (or towards 
filling up the dent), creating a net current towards the 



region of larger height, i.e. uphill as shown by the thick 
white arrow in Fig. (a) . 



a) lower terrace ~i ll » lower terrace c) lower terrace 
higher ten-ace upper terrace ^ higher terrace 



FIG. 15. Schematic illustration of the instability caused 
by the line tension in the step-edge. View from above, the 
thick line is the boundary between two terraces whose differ- 
ence in elevation is one monoatomic layer, a) shows a contin- 
uum picture while b) and c) represent the same on a lattice 
for different step-boundary orientations. 

Note that although this current stabilizes the step- 
edge, it destabilizes the 2+1 dimensional surface h(x,y) 
at long wavelengths, creating mounds on the substrate 
induced by the uphill current. The local current density 
Ji along the step-line is proportional to the local gradient 
along the edge of the chemical potential /x, i.e., one can 
write Ji = — ^V e /i, where v is the mobility of the atoms 
along the step edge [^2| . Assuming local thermodynamic 
equilibrium, the chemical potential may be expressed as 
a functional of a free energy which in turn is written 
as that of an clastic line of line-tension a. The result 
is that /i is proportional to the local curvature k of the 
\. k is the local curvature on the line 



line, /i = — (77 

h(x, y) = const, and thus it can be expressed in terms of 
the local derivatives of h: 



Ji = -o-vV \Vh\- 
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(9) 



Next, we calculate the uphill current due to the cur- 
rent density (in Eq. (||)) created by an arbitrary shaped 
hump created at the edge of a single straight step, and 
then specialize the expression for a simple, cosine shaped 
hump. We assume that one can choose an (x, y) coor- 
dinate system such that the equation for the step-profile 
h(x,y) = const in this system is simply expressed by a 
single valued function: y = y(x), see Fig. [l^. In this 
simple setup the edge is parametrized by x, and thus the 
step-edge current in Eq. (0) becomes: 



Ji = va - —^tb ( e * + V e y 



(10) 



where y^\ y" and y' are the third, second and first order 
derivatives of the profile, and e x and e y are the unit 
vectors along the x and y axes, respectively. The total 
mass transported per unit time by curvature gradients 
between points A and B on the step-edge is defined as 
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the line-integral of the edge current density: 

Iab = I ds Ji 

J A 



(11) 



where ds is the integration element along the curve of the 
step-edge. In the following we illustrate that this current 
points toward the base of the higher step. 
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FIG. 16. A straight step-edge with a hump. The step-edge 
current density generated by the local curvature gradient 
along the hump creates a net uphill current. 

Let us choose a step-profile given by y(x) = a cos (bx) 
between the points xa = —it/b and ig = n/b. After 
performing the integrals in Eq. (Ill]), one obtains: 



l AB 



Avab 
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I + H2 + 6z z 



(l + z) 2 
1 + 3z 



l + z 



(12) 



where z = a 2 b 2 is a dimensionless number, and K(fc) 
and E(fc) are complete elliptic integrals of first and sec- 
ond kinds respectively. Since l + llz + 6z 2 > l + 4z + 3z 2 , 
and E(iy/z) > K(iy/z), for all z > 0, the mass-current 
has a negative y component (the x component is zero 
because of symmetry), i.e., it points towards the base 
of the higher step, uphill. In the limit of small humps, 
z -C 1 (small slopes approximation), K(jy / z) ~ §), 
E(iy/z) ~ \ (1 + §), so the current becomes proportional 
to the square of the hump's height, Iab ~ —ixvob 3 o? 
and for peaked humps (z ^> 1), K(i^/z) ~ In (iy/z)/y/z, 
E(i-y/5) ~ yfz the current depends linearly on the height 
of the hump, Iab ~ \uvb 2 a. The same analysis can be 
repeated when there is a dent in the step profile, with 
similar conclusions. Certainly, the currents expressed 
above are instantaneous local currents. If one is inter- 
ested in the evolution of the step profile, one can employ 
simple geometrical considerations (see p3,E3|) to write: 



dy 
dt 



(13) 



where y = y(x,t), v n is the normal step-edge veloc- 
ity and the prime denotes derivative with respect to x. 



Assuming that the mass transport is solely due to the 
step-edge current, one has v n ds = —(dJi/dx)dx, with 
ds = dxyjl + (y') 2 . Thus Eq. @ takes the form: 



dy 
dt 



dJj_ 

dx 



y 



(4) 



10y'y"yW 
~WW ~ [1 + (y') 2 } 3 
3[l~5(2/') 2 ](y") 3 



[i + (y r 



214 



(14) 



In the limit of small slopes, y' << 1, y" << 1, y^ << 1, 
?/ 4 ) << 1, and so from Eq. ( |l4| ) it follows that: dy/dt = 
— vay( A \ after keeping the leading term only. This is the 
well known Mullins equation for the dynamic relaxation 
of a step-edge, related to the fourth order linear growth 
equation, Eq. (||) with A22 = 0, followed by the LC model 
in our growth simulations. 

The continuum description presented here com- 
pcllingly demonstrates the possibility of a destabilizing 
uphill current arising entirely from step edge diffusion — 
an SED instability without any ES barriers, as discussed 
111 Refs. |,j§ recently. 



B. A discrete description 

In reality the deposition process (and our simulation of 
limited mobility models) takes place on the discrete and 
atomistic crystalline lattice, which introduces an orienta- 
tion dependence of the line tension a. One would natu- 
rally expect that the high-symmetry, in-plane crystalline 
directions will have the largest a, these being the most 
stable. However, there is a hierarchy even among these 
high-symmetry orientations, as is illustrated in Figs. [l5| 
(b) and [l5| (c). Fig. |l5| (b) shows a step-edge aligned 
along a diagonal of the square lattice and in Fig. [15] (c) 
the step-edge is oriented along one of the main axes of the 
lattice. While the step is stable along the diagonal, it is 
not as stable along the main axis, since in this latter case 
(of Fig. [if] (c)) in order for an atom to reach a higher 
coordination site along the line, it has to detach itself 
from the step-edge (to become an adatom) first, which is 
energetically less favorable. 

Let us now analyze in more detail the effects of SED 
in the case of the two most important MBE-motivated 
limited mobility models, WV and DT models. The dif- 
fusion rules for the mobile atoms for these two models 
are only slightly different: in the case of WV model par- 
ticles seek to maximize their coordination while in the 
DT model they only seek to increase it. Nevertheless 
the two models present completely different behaviors. 
As we have seen in Subsection II. C, the NRT version of 
WV model produces mounded structures, however this 
does not happen for the DT model. In the following we 
explain this unexpected difference between WV and DT 
models based on local SED currents. Given the rather 
"minor" differences between DT and WV rules (and their 
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almost identical growth morphologies in 1+1 dimensions) 
it is remarkable how different their 2+1 dimensional NRT 
growth morphologies are (WV has spectacular mounded 
morphology as seen in Fig. [| and DT has essentially sta- 
tistically scale invariant kinetically rough growth as in 
Fig. @). 



a) "WV" 



lower terrace 



(110) 



r 



rl 




j 

i) 

- do not compensate 

- "net current is uphill" 

do not compensate 

, compensate 

- J upper terrace 

compensate 

b) "WV" lower terrace 

compensate are not compensated ( 1 00) 




upper terrace . .... 

r * net uphill current 



FIG. 17. Local currents along step-edges of orientations 
(110) and (100) for the WV rules. 

Providing an understanding of this difference (between 
DT and WV morphologies) in terms of the SED instabil- 
ity (WV has it and DT does not) is one of the important 
accomplishments of our work. 

Figs, [l?] and [l^ show two high-symmetry oriented 
steps, (110) and (100) respectively, with a small dent 
in it. The local particle current shown by arrows are cal- 
culated using the rules of WV and DT models. One can 
immediately conclude that in the case of these configu- 



rations, the net current is jWX) = and j 



wv 

'(100) 



= 2 



both directed uphill, inn) = ~~ v^/2 directed downhill 
and in oo) =1/3 directed uphill. (The current is consid- 
ered to be 1 across the boundary between nearest neigh- 
bor sites X and X', if a particle once at position X will go 
in one step with probability 1 to position X'. In case of 
ties, the probability fraction corresponding to the rules 
is taken). 
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FIG. 18. Local currents along step-edges of orientations 
(110) and (100) for the DT rules. 

As seen from the configurations above, there is a clear 
net uphill local current in the case of WV model in both 
cases of high-symmetry orientations. According to the 
theory of SED, one should see a mounded morphology. 
Nevertheless, when we look at the morphology of Fig. 
^| (a), one only sees a rough surface not obviously a 
mounded landscape in the n r = 1 WV growth. The rea- 
son for this discrepancy lies at the heart of the SED insta- 
bility and the conditions (and constraints) for its observ- 
ability: having destabilizing step-edge diffusion currents 
is not sufficient for the creation of mounds. There is an- 
other very important necessary ingredient (which was not 
analyzed in Refs. |||||) that is needed for the formation 
of quasi-regular mounded growth patterns in addition to 
the destabilizing SED currents: suppression of deposition 
and nucleation noise. In order to have a consistent uphill 
current build-up (to create mounds) , we need a sustained 
stability of the step-edge, in other words, if we look at 
the step-edge as a 1+1 dimensional 'surface' in the x — y 
plane (the upper terrace being the 'substrate') we need 
to have more or less a solid-on-solid type of step-edge, 
with perhaps some small overhangs eventually, but no 
large or frequent ones. The only effect which can disrupt 
a stable step-edge is the nucleation of new islands on the 
terraces near the edge. Let us imagine that we start with 
a straight step, aligned along the most stable direction. If 
island nucleation takes place easily, then the islands that 
nucleate close to the stable edge will quickly grow into the 
step, disrupting the step's contour and identity. If the re- 
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laxation time of the disrupted contour is larger than the 
nucleation time, small newly grown islands will disrupt 
the step even further until it ceases to exist. In this case 
the SED instability will not manifest in the mound for- 
mation because the necessary condition for mound for- 
mation is not satisfied globally although there may be 
local uphill current. Thus, noise may prevent the SED 
instability from manifesting a global mounded instability. 

The importance of noise suppression for the SED in- 
stability is strikingly apparent from comparing Figs. |B| 
(a) and | (b) for the WV model. The NRT version has 
mounded morphology and the n r = 1 version does not. 
However for the DT model, one finds no such regular 
mounded structures even when noise and nucleation is 
strongly supressed (Fig. ||). The reason for this differ- 
ence is that in the DT model the SED current is essen- 
tially downhill (and thus is stabilizing). Depending on 
the edge configuration (for example the (100) edge in 
Fig. |l8| ) one may find a small uphill current in the DT 
model, however it is statistically not significant to desta- 
bilize the whole surface. There will be a weak mounding 
tendency at early times, when the V 4 /i diffusion term is 
relevant and when the height-height correlation function 
shows some oscillations (see 0), indicating the presence 
of this term, but then it will be quickly taken over by the 
stabilizing downhill current and nucleation events. Thus, 
the DT model, at best, will exhibit some weak and irreg- 
ular epitaxial mounding Q , but not a regular patterned 
mounded morphology. 

Another important point is that mounding due to the 
SED instability does not require any growth nonlinear- 
ity to be present. This is obvious from the fact that the 
LC model, which is strictly linear by construction, shows 
spetacular epitaxial mounding induced by the SED in- 
stability. In the following subsection we analyze in more 
details the effects of noise reduction on the SED instabil- 
ity. 

C. Noise reduction method: a numerical way to 
control the step edge diffusion instability 

As emphasized in the previous subsection, a necessary 
ingredient for creating mounds by SED is the suppres- 
sion of noise. In the following we restrict our analysis to 
discrete limited mobility models, and frequently use for 
comparison WV and DT models. Unfortunately, rigorous 
analytical calculations are practically impossible for dis- 
crete and nonlinear dynamical growth models of DT and 
WV types. However, we make an attempt by introducing 
the proper quantities to pinpoint the effects of noise re- 
duction on these two models and actually show the major 
difference between the two models in one and two sub- 
strate dimensions in the context of SED instability. First, 
we introduce the notion of conditional occupation rates. 
The conditional occupation rate /; associated with site 



i in a particular height-configuration is the net gain of 
transition rates coming from the contribution of all the 
neighbors within a distance of I sites on the surface and 
from the growth direction (from the deposition beam). 
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FIG. 19. Conditional occupation probabilities for two par- 
ticular configurations in 1 + 1 dimensions for DT and WV 
models. The fi-s that differ for the DT from WV are set in 
bold italics. Here 1 = 1. 

The fi-s can easily be calculated once the diffusion 
rules are known. If Wj_>j' is the transition rate for an 
atom landed on site i to be incorporated at site i' (which 
is a nearest neighbor when 1 = 1), then: 

{<}l 

f t =1+^(Wm-Ww) (15) 

v 

where the unity in front of the sum represents the con- 
tribution from the beam and {i}i means summing the 
contributions over all the neighbors of site i within a dis- 
tance I. The fi-s are obviously conditioned to the event 
that there is an atom landed on site i or on the neighbors 
i' . It is important to emphasize that these quantities in 
Eq. (|TJ) are introduced for the no-desorption case only: 
once an atom is incorporated in the surface it stays that 
way. If one wishes to add desorption to the model, then 
the evaporation rates have to be included in Eq. (|l5| ) as 
well (the rates in Eq. (|l^) are all diffusion rates ). 

Let us now illustrate in light of the two models WV and 
DT, the conditional occupation rates, in 1+1 and 2+1 di- 
mensions. Figs. |l!^ and ^ show the fi-s for particular 
configurations for one and two dimensional substrates re- 
spectively. The fractional transition rates mean that the 
motion of the atom is probabilistic: it will choose with 
equal probability among the available identical states. 
This is in fact the source of the diffusion noise r\ Cl aris- 
ing from the stochastic atomistic hopping process, which 
should be added to the right hand side of Eq. (|l]), how- 
ever, it is an irrelevant contribution in the long wave- 
length scaling in the renormalization group sense. The 
conditional occupation rates fi are small rational num- 
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bers from the interval: 



(d) 



(16) 



where Z\ is the maximum /th order coordination number 
on the d-dimensional substrate. 
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FIG. 20. Conditional occupation probabilities for two par- 
ticular configurations in 2 + 1 dimensions for DT and WV 
models. The fi-s that differ for the DT from WV are set in 
bold italics. Here 1 = 1. 

For example, on a cubic lattice , — 2d. For a hxed 
set of deposition and diffusion rules the number of possi- 
ble occupation rates is hxed, and for low Z-values it is only 
a few. This number is determined by how many differ- 
ent values can be obtained from all the possible relative 
height combinations of all sites within a ci-dimensional 
ball of radius 21 + 1 (the gain coming from the neighbor 
at Z-away may depend on the configuration of radius I 
around that neighbor) . If I = 1 then in 1 + 1 dimensions 
this means that a chain of maximum 7 sites has to be 
considered and in 2 + 1 dimensions it is a square domain 
of maximum 25 sites to exhaust all the possibilities (the 



domain in this latter case is given by all sites that 
obey \i\ + \j\ < 3). For example, the one dimensional 
WV model will have conditional occupation rates from 
the set {0, 1/2, 1, 3/2, 2, 5/2, 3}. For the one dimensional 
DT, the set is a bit restricted: {0, 1, 3/2, 2,5/2, 3}, it does 
not allow for / = 1/2. 

Fig. |2l] shows the plot for the occurrence probabilities 
P(f) versus the conditional occupation rates for both the 
WV and DT models (n r — 1). The occurrence proba- 
bility of rate / in one dimension is calculated as follows: 
take all the possible height configurations that are distin- 
guished by the deposition and diffusion rules, (e.g. WV 
and DT) of the 2(2Z + 1) + 1 sites, monitor the / that 
belongs to the middle site and then measure the rate of 
occurrence of each value for /. 




FIG. 21. Conditional occupation rates and their occurence 
probability p(f) for the WV and DT models. Here 1 = 1. 

It is important to count these occurrences relative to 
the distinguishable number of configurations. Both WV 
and DT models have their rules for the motion of the 
landed atom formulated in terms of coordination num- 
bers. In the WV model, the landed atom will move to the 
site with the maximum coordination whereas in DT the 
atom will move to increase its coordination. Thus only 
those configurations will be distinguished by the rules 
that have different coordination numbers in at least one 
of the sites. A given set of coordination numbers thus 
determines a class of configurations. For example in 1+1 
dimensions both the WV and DT models on 7 sites will 
distinguish a total of 3 6 = 729 classes, since the depo- 
sition and diffusion rules for these two models do not 
involve actual height differences, just the fact that a site 
is at higher, the same, or lower level than its neighbors. 
In 2 + 1 dimensions the exact enumeration of these classes 
is highly non-trivial, and it constitutes the subject of a 
future publication. 

One can note the following property of conditional oc- 
cupation rates: in every domain of sites on the substrate 
that has a zero inflow and zero outflow current across 
its boundaries ('island'), the sum of the rates within the 
island adds up to the total number of sites within the 
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island: 



E 

{^island 



fi = nr. of sites within the island. 



(17) 



This is a consequence of conservation of probability. The 
span of the island may only change if a new particle is 
deposited within a distance 21 + 1 from its boundary 
outward from the island. This is valid in any dimension. 

The larger the fi the more probable it is that site i 
will receive an atom (it can come from more directions 
than for one with a lower /) . We believe that when using 
the NRT, i.e., n r > 1, the sites with higher / will have a 
higher chance to be actually deposited onto. To illustrate 
this, we consider the simplest setup: starting from a fiat 
surface with a few irregularities, we perform the counting 
process and see with what probability the various sites 
will actually be deposited onto when the corresponding 
counters reach the threshold value n r . 
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FIG. 22. Particular counter values with different noise 
reduction parameters a) n r = 4, and b) n r = 12. c) shows 
the probability that site i reaches the treshold for the first 
time, d) Plot of probability for nucleation on the terrace 
(of L — 3 sites) and deposition at the dent site with / = 3. 
The straight line is an exponential fit with slope —0.285 and 
intercept 8. 

Unfortunately rigorous analytic derivations are diffi- 
cult even in the simplest case because of the fact that 
surface growth is driven by the first-passage time proba- 
bilities of the counters (or waiting times). We therefore 
follow numerically the evolution of counters on a sim- 
ple substrate which is flat everywhere except at one site 
where the height is lower by —1 (see the surface profile 
just inder Fig. ^2] (a)). According to both rules (WV 
and DT) the occupation rate is 1 except for the lower 



site where it is 3 and its two nearest neighbors where 
it is 0. Fig. ^ (c) shows the probability that the first 
real deposition event occurs at site i for n r = 8, L = 40. 
Obviously all sites on the flat region will have the same 
chance of P(f, n r ) = P(l, n r ) to reach the threshold for 
the first time, while the site with / = 3 will be the most 
probable to reach it, and for its two nearest neighbors it 
is entirely unlikely: P(0,n r ) = 0. In case of Fig. ^2] (c), 
P(l,8) = 0.0077 and P(3,8) = 0.71. This means that in 
spite of the uniform randomness of the shot noise, 71% 
of the time the actual deposition will occur on that single 
site with higher coordination and only 21% of the time 
we will actually be able to land an atom somewhere on 
the whole flat island of 37 sites. This clearly shows that 
nucleation of new islands on flat regions is strongly sup- 
pressed. In Fig. [22] (d) we represent the two probabilities 
P(3,n r ) and P(l,n r ) as functions of n r (on a log-linear 
scale). For higher noise reduction parameter values we 
are practically forced almost each time to fill in a site 
with the highest occupation rate, and with almost zero 
probability to land an atom on the flat terrace. Looking 
back on Figs [l?] (a) and [H] (b) for the 2+1 WV model, it 
becomes clear that the main particle current will be along 
the step edge (the shaded squares on the figures) and it 
will be greatly enhanced by NRT. At the same time the 
rate of island nucleation is drastically reduced, so the re- 
laxation of the step-edge becomes rather fast compared to 
the creation of new islands. Since the step-edge dynamics 
becomes fast compared with the nucleation dynamics the 
mounds that are formed will grow into each other pre- 
serving their pyramidal shapes, exhibiting the observed 
mound-coarsening behavior. Models in 2+1 dimensions 
can be related to the conditional occupation rates shown 
in Fig. ^0|. The hole in the middle of the top island, and 
also the dents in the step-edges have much higher rates 
to be occupied in the WV model than in the DT model. 
Thus the step-edge currents are a lot stronger for the WV 
model than for the DT model. In addition the middle site 
at the base of the upper island has a lower occupation 
rate for the WV (/ = 2) than for the DT (/ = 5/2 = 2.5) 
model. This means that in the case of WV it is harder to 
disrupt a straight step-edge than for the DT model. We 
believe that the spectacular difference between the DT 
and WV noise reduced morphologies discovered in this 
work arise from the conditional occupation differences in 
the two models, which lead to an SED instability in the 
WV, but not in the DT model. While we have provided 
a rather compelling qualitative picture in this subsection 
explaining why the SED instability produces spectacular 
mounded morphologies in the NRT version (n r ^ 1) of 
the WV model, but neither in the DT model (n r = 1 or 
n r =/= 1 version) nor in the n r = 1 WV model, we have 
not been able to come up with an analytical theory for 
the mound formation, which remains an important open 
problem for the future. 
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IV. PROPERTIES OF THE MOUNDED 
MORPHOLOGIES 

In this section, we discuss the dynamical scaling prop- 
erties of the mounded morphologies presented in Section 
II. A conventional way to study mounding is to 

look at the height-height correlation function 



H(r) = (/i(x)/i(x + r)) x , 



(18) 



where r — |r| is the distance between two sites on the 
substrate, and the (...) denotes a substrate averaging. 



order (LC model), the nonlinear fourth order, and the lin- 
ear sixth order models. Note that noisy oscillations are 
already detectable in the original models with n r = 1, 
but the oscillations here are weak and irregular Q. The 
oscillations become spectacularly enhanced and rather 
evident when the noise reduction technique is used. This 
corroborates the fact that NRT suppresess the noise, but 
does not change the growth universality class. 

In Fig. [m] we present the non-oscillatory H{r) corre- 
sponding to the non-mounding interfaces of F, DT, and 
the nonlinear fourth order equation with infinite order 
nonlinearities (Eq. (R)) models . 
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FIG. 23. The height-height correlation functions H(r) 
corresponding to the mounded morphologies presented in Fig. 
H |, |, and O. From top to bottom : the WV model at 10 4 
ML, the LC model at 10 3 ML, the nonlinear fourth order at 
10 3 ML, and the linear sixth order at 10 6 ML. 

The calculated H(r) oscillates as a function of r in sys- 
tems with mounded patterns. Fig. shows the oscillat- 
ing correlation function from the WV, the linear fourth 
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FIG. 24. The height-height correlation functions HM 
corresponding to the morphologies presented in Fig. ^, ^, 
and |lo[ From top to bottom : the F model at 10 6 ML, the 
DT model at 10 6 ML, and the nonlinear fourth order with 
higher nonlinearities at 10 3 ML. 

The morphologies from the original version (n r = 1) of 
these models are kinetically rough. But when the noise 
is suppressed with n r > 1, no specific oscillatory pattern 
emerges although some weak and irregular H(r) oscil- 
lating could sometimes be detected in the DT model Q 
(particularly for n r = 1), most likely caused by the pres- 
ence of the V 4 h Mullins term in the DT growth equation. 
As a matter of fact, these rough surfaces (F, DT, etc.) 
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merely become "smoother" and flatter for n r > 1. 
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FIG. 25. The average mound radius as a function of time 
from the WV model with n r = 5 from a substrate of size 
100 x 100. The average mound slope of the same system is 
shown in the inset 

Other important characteristics of mounded patterns 
P,pT| are the average mound radius and the average 
mound slope. Conventionally, the distance of the first 
zero-crossing of the correlation function H (r) is taken to 
be the average mound radius R, and [H(r = 0)] 1 / 2 is the 
average mound height. 
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FIG. 26. The average mound radius as a function of time 
from the LC model with n r — 5 from a substrate of size 
100 x 100. The average mound slope of the same system is 
shown in the inset 

The average mound slope (M) is naturally the ratio of 
average mound height and average mound radius. The 
scaling relation for R can be written as 

R[t) ~ i" (19) 

where n is the coarsening exponent. The average mound 



slope M increases in time as 

M(t) ~ i A 



(20) 



where A is the steepening exponent. In Fig. [2^ we show 
the evolution of the average mound radius in the WV 
model with n r = 5. The average mound radius is ob- 
viously increasing in time, showing that the coarsening 
process is dominant during the growth period with the 
coarsening exponent n = 0.24. The inset shows the av- 
erage mound slope which does not change much in time 
(A~0). 
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FIG. 27. The average mound radius as a function of time 
from the nonlinear fourth order model with n r = 5 from a 
substrate of size 100 x 100. The average mound slope of the 
same system is shown in the inset 

This shows that the mounds in the WV model have 
constant magic 'slopes' or equivalcntly, the WV mounded 
growth exhibits slope selection. The value of this selected 
slope in the WV model is approximately tan 9 w 0.1—0.2 
or 8 w 6 — 12° and it does not seem to depend strongly 
on the value of n r . Similar behavior is also found in the 
SED-induced mounded morphologies in the LC model 
with n r = 5 (Fig. |2^) and the nonlinear fourth order 
model with n r = 5 (Fig. |?]). We note that the scaling 
properties of the ES barrier-induced mounded morpholo- 
gies are quite distinct from the SED-induced mounding 
instability. Asymptotically, there is no slope selection in 
the mound formation caused by an ES barrier f2~l]| . Af- 
ter an initial crossover period, the coarsening drastically 
slows down with n — > and the steepening exponent 
becomes exceptionally large with A — + 0.5. Thus SED 
induced mounding has slope selection and ES barrier in- 
duced mounding typically does not (asymptotically). 



A. Comparison with experiments 

In this subsection we discuss experimental epitaxial 
mound formation in light of our findings about the lim- 
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ited mobility growth models. First and foremost we cau- 
tion against the wisdom of automatically assuming the 
existence of an ES barrier whenever a mounded morphol- 
ogy is seen. As we have argued in this paper, the ES 
barrier is most certainly a sufficient condition to create 
mounded growth morphologies (nominally without any 
slope selection, unless some additional mechanisms, not 
intrinsic to ES barrier, is invoked), but by all means, it 
is not the only possible cause or a necessary condition. 
When analyzing experimental mounding data, therefore, 
an ES barrier should not be automatically assumed and 
more than one possibilities should be critically consid- 
ered. We hope to make our point by using a recent ex- 
periment as an example . It is instructive to carry out 
a comparison between the observed mound morphology 
in Cu(100) growth pi[ and the mound morphology from 
the LC (the simple linear fourth order) model shown in 
this paper. 
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that the simple LC model without ES barrier can produce 
mounded morphology that is almost exactly the same as 
the Cu(100) mounded growth and the possibility that the 
mound morphology in Cu(100) may involve mechanisms 
other than the ES barrier should be explored. Clearly the 
SED mechanism as a possible source for slope selected 
mounding in Cu(100) growth should be taken seriously. 

Finally, we recall a "puzzling" observation from Ref. 
p4j , namely, that the average mound slope calcu- 
lated from the ratio between the average mound height 
([H(r = 0)] 1 / 2 ) and the average mound radius (first zero 
crossing of H(r) correlation function) is much smaller 
than the typical mound slope measured from the mor- 
phology directly. In the Cu(100) case J24|, the average 
selected slope is 2.4° and the typical mound slope is 5.6°. 
In our LC study, the average slope is approximately 6° 
while the typical slope is in the range between 10° to 15°. 
This is so because [H(r — 0)] 1 / 2 (taken as the average 
mound height) does not equal the typical mound height. 
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FIG. 28. The contour plot of the "equal priority" 
morphology with I = 2 and n r = 8 at 10 4 ML. 
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In Fig. E8[ we show a contour plot of our mounded 
growth front simulated from the simple LC model with 
diffusion length I = 2 and noise reduction factor n r = 8. 
The plot looks very similar (in fact, almost identical) to 
the STM images of the mound morphology in reversible 



homoepitaxy on Cu(100) in Ref. 24|; both systems show 
mounds with square bases and the mounds are lined up 
in a somewhat regular pattern. Not only are the two 
systems strikingly comparable visually, but the quanti- 
tative dynamical scaling properties for the two cases are 
also approximately the same. In the experiment [Q, the 
mounds have a magic slope of 5.6 ± 1.3° and the coars- 
ening exponent is given as n = 0.23 ± 0.01. In our LC 
simulation, n w 0.25 and the slope selection process is 
clearly observed (see Fig. |2^) with typical mound slope 
being around 10°. The selected slopes in the two sys- 
tems are not equal, but this is not particularly crucial. 
The LC model is only a simple dynamic growth model 
which in general should not be compared quantitatively 
with real MBE growth. We are not making a claim that 
there is absolutely no ES barrier effect in the Cu(100) 
epitaxial growth. What we want to point out here is 



FIG. 29. A simple shematic setup to show the difference 
in the definitions for the average mound slope. 

Let us consider the simplest "mounded" pattern which 
mimics the height variation in the morphology along 
a cut perpendicular to the base of the pyramids (see 
the schematic Fig. ^). In the simplest situation this 
is described by the periodic, piecewise linear function: 
h(x) = -^(x- AkR) + A, if x e (AkR,Ak + 2R), and 
h(x) = ±[x-4.kR-2R)-A, if x S (AkR + 2R, Ak + AR), 
k = 0,1,2,..., where the meaning of A and R is obvi- 
ous from Fig. Thus, one has for this simple shape 
[H(r = 0)] 1 / 2 = (A/R)/VH = 0.57.. s, where s = A/R 
is the typical mound slope, illustrating the difference be- 
tween the two definitions. In addition we also have a 
small discrepancy introduced by the fact that slope mea- 
sured from the ratio of two averages is different from 
the average of the ratios. Thus, our simple LC model 
mounding explains another puzzling aspect of experi- 
mental mounding not easily explained via the ES bar- 
rier mechanism. We refrain from carrying out further 
comparison with other experimentally reported epitaxial 
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mounding in the literature because the main point we are 
making is qualitative: SED may be playing an important 
role in observed epitaxial mounding. 

B. Global particle diffusion currents 

In the conserved growth models where the continuity 
equation, Eq. (||), applies, global particle diffusion cur- 
rents can be measured directly from growth on tilted 
substrates |25[ . If the measured current is a downhill 
current (negative) then it corresponds to a stable sur- 
face that belongs asymptotically to the EW universality 
class, and an uphill particle current will point to an unsta- 
ble surface (i.e. mounded morphology) p5j . Our parti- 
cle current results, however, show that this concept may 
be misleading for growth in d=2+l dimensional physi- 
cal surfaces since the global current may show no simple 
systematic behavior as a function of tilt. In general, sub- 
strate current in 2+1 dimensional growth may be positive 
or negative depending on tilt direction. In fact, we be- 
lieve that the current measurement on tilted substrates is 
a useful technique for determining the universality class 
of growth models only in 1+1 dimensions. In 2+1 di- 
mensions, where the surface current could be stabilizing 
in one direction (e.g. 100) and destabilizing in another 
(e.g. Ill), the technique fails. This was not appreci- 
ated in Ref. |^5| and has not been earlier discussed in the 
literature. 

An important example is the particle current in the 
WV model in various dimensions. It was reported |25| ] 
that the WV model has downhill current in both one 
and two dimensions. This implies that the WV model 
belongs to the EW universality class, which was the con- 
clusion in Ref. ^5|. However, our Fig. |j] shows that the 
two dimensional version of the WV model has mounded 
morphology induced by an SED instability. This seems 
confusing as mound morphology should have uphill cur- 
rent. We repeated the particle current study on the two 
dimensional WV model and found the same downhill cur- 
rent as in the literature pq| , both in the original WV 
model with n r — 1 and in the n r = 5 WV model. But all 
the global particle current measurements in the literature 
were obtained (2^j2(| by tilting the substrate in the (100) 
direction! We have also studied the diffusion current by 
tilting the substrate along the diagonal, i.e. in the (111) 
direction , and found a completely different result. The 
WV current from the (111) tilted substrate is uphill (pos- 
itive). In the n r = 1 model, the (111) current starts out 
as uphill and then crosses over to a negative value, i.e., 
downhill with crossover time increasing with increasing 
substrate size L. This is already a weak indication that 
there are some instabilities present. In the n r > 1 WV 
model, however, we measure strong (111) uphill current 
which stays positive up to 10 6 ML which is our longest 
simulation time, confirming the suspicion that in fact the 



WV model in two dimensions has the SED mounding in- 
stability induced by (111) uphill current, and a simple 
measurement of (100) downhill current gives a qualita- 
tively incorrect result. 

Our result shows that the global particle current 
strongly depends on how the substrate is tilted. This tilt 
direction (e.g. 100 versus 111) dependence may be very 
complicated in systems with two or higher dimensional 
substrates since surface current is a vectorial quantity 
in 2 (or higher) substrate dimensions. The particle cur- 
rent from one particular tilt may be highly misleading 
with respect to the overall behavior of the system. The 
downhill (100) current in the WV model had incorrectly 
suggested that the model belongs to the EW universality 
class and should exhibit relatively flat morphology in two 
dimensions. However, by properly suppressing the noise 
our study shows that the two dimensional WV model 
actually becomes unstable due to a strong SED driven 
uphill (111) current. As a matter of fact, there has been 
at least one prior report of mounding in the WV model 
in unphysical higher dimensions (three and four dimen- 
sional substrates) pTj]. In higher substrate dimensions 
the noise effects are intrinsically a lot weaker than in one 
or two substrate dimensions, because diffusion in lower 
dimensional substrates is topologicaly more constrained, 
and therefore the SED instability is easier to observe 
(even without noise reduction) in higher dimensions. We 
have explicitly verified through direct numerical simula- 
tions that higher (3+1 and 4+1) dimensional WV growth 
has spectacular SED instability-induced mound forma- 
tion even in the absence of any noise reduction, thus ver- 
ifying the old and unexplained results of Ref. ^7|. The 
higher dimensional DT growth is smooth and flat, and 
exhibits no such mounding. This arises from the much 
stronger edge diffusion effects in higher dimensions which 
effectively introduces very strong intrinsic suppression of 
noise. Our work therefore resolves an old puzzle, ex- 
plaining why the higher dimensional simulations |27| of 
the WV model produced mounded morphologies. 

V. CONCLUSIONS 

We have shown that one has to carefully account for 
subtle topological-kinetic instabilities when trying to un- 
derstand mounding in the surface growth process. These 
instabilities are simply the result of the fact that in higher 
(higher than 1+1) dimensions a local current may act as 
a stabilizing one in a certain projection on the substrate 
and as a destabilizing one in another projection. The lo- 
cal anisotropy of current effects averaged over the whole 
substrate may generate mounded growth even in situa- 
tions (e.g. WV model) where mounding is completely 
unanticipated. The situations where SED instabilities 
may be present are rather subtle to analyze and each 
case has to be investigated rather carefully using the 
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noise reduction technique. The fact that SED mound- 
ing instability exists in the noise reduced WV, LC, and 
the nonlinear fourth order continuum model, but not in 
F, DT, and the infinite order nonlinear model shows that 
the kinetic-topological nature of SED mounding may not 
be obvious at all in a given situation. We note here that 
the DT model also exhibits very weak mound-like mor- 
phologies f|. However, the weak mound- like structures 
in the DT model arise from the fact that the roughness 
exponent of the model is exceptionally large (arising from 
the presence of the V 4 /i term in the growth equation). 
The mounds (and the oscillation in the height correla- 
tion H(r)) in the 2+1 DT model are random and irreg- 
ular while the SED-induced mound formations (in e.g. 
noise reduced WV model) exhibit regular patterns with 
strong periodic oscillations in H(r). 

One of our important findings is that SED gcncrically 
leads to slope-selected mounded morphology whereas 
the ES instability generically leads to slope steepening 
mounding. We believe that many of the slope selected 
mounding morphologies observed in MBE experiments 
(with typical selected slopes being small ~ 10° or less) 
may actually arise from the kinetic-topological mound- 
ing we find generically in the LC or the WV model in 
this paper. Further experimental and theoretical work is 
required to settle the competing roles of various instabil- 
ities contributing to epitaxial mounding in real surface 
growth. 
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